Boundary height fields in the Abelian sandpile model 
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■ Abstract 

We study the abelian sandpile model on the upper half plane, and reconsider the correlations of the 
^ , four height variables lying on the boundary. For more convenience, we carry out the analysis in the 

I dissipative (massive) extension of the model and identify the boundary scaling fields corresponding to 

' the four heights. We find that they all can be accounted for by the massive pertubation of a c = —2 

0^ ' logarithmic conformal field theory. 

o ■ 
^ : 

^ ; 1 Introduction 

^ , The description of equilibrium critical phenomena has been one of the greatest success of two-dimensional 

Q I conformal theories in the past twenty years PP. More recently, attention has focused on new types of 

' obscrvables in otherwise well-known models, and also on new kinds of critical systems. In both cases, 

^ ■ non local features often play an important role, because either the observables one is interested in are 

. themselves non local in terms of the natural microscopic variables, or else because the statistical model 

possesses intrinsic non local properties. These studies often lead to a description in terms of conformal 
theories with peculiar properties. A class of systems with such properties is provided by sandpile models. 
Some of them, and in particular the one we consider here, are believed to have a faithful description in 
terms of logarithmic conformal field theories. 

Our motivation to study these models is two-fold. First, one wishes to see to what extent they lend 
themselves to a conformal field theoretic approach, and if the adequacy of the conformal description 
is as good as for the equilibrium systems. Second, the logarithmic theories have been developed for 
themselves, but are complex and some of their aspects are not fully understood yet. It should therefore 
be profitable to have concrete realizations in order to have a better understanding of the most peculiar 
features. 

The sandpile model we consider here is the isotropic Abelian sandpile model (ASM), as originally 
defined in [2]. The most natural variables to consider in a conformal context are the four height variables. 
In the bulk, correlations of height 1 variables can be handled by local calculations [3], but height 2, 3 
and 4 variables are much more complicated, and only their one-site probabilities are known For 
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sites on a boundary, closed or open, Ivashkevich showed, by using suitable identities, that the non local 
configurations needed to handle the heights bigger than 1 could be reduced to local computations 
He was then able to compute the two-site joint probabilities of all height variables. He found that all 
correlations decay like r~^, and infered, wrongly as we will see, that all boundary height variables scale 
to the same conformal field. 

Our aim in this paper is to revisit this problem and to identify the fields corresponding to the 
boundary height variables, but in a dissipative extension of the model, known to be described by a 
massive perturbation of the c = — 2 conformal theory [Hj. The main advantage for doing this is that it 
allows for an unambiguous identification of the fields from a few 2-point correlators, off, and therefore 
also at, criticality. The so-obtained identification can then be checked from other 2-point and from 
3-point functions. In contrast, when one considers the critical, non dissipative model, the 2-point 
functions yield ambiguous field identifications, which can only be fixed by using 3-point functions, and 
then checked from higher correlators. 

The article is organized as follows. The next section defines the model and sets our conventions. 
Section 3 deals with the boundary unit height variable and what we call supercritical height variables. 
They are much easier than the other ones, and illustrate the way the identification with concrete fields 
is obtained. The identification of the height 1 and supercritical height fields also facilitates that of the 
other heights. 

Sections 4 and 5 form the hardcore of the article. In Section 4, we explain our prescription — a 
two-step burning algorithm — to associate recurrent configurations of the sandpile to spanning trees, 
from which a clear characterization of the local height constraints follows. We use it to compute 2-site 
probabilities for having a height 1 or a supercritical height at one site, and any other height at the other 
site, from which we deduce the field identification of all height variables. Section 5 checks these results 
by computing explicitely all 2-site height probabilities in the massive extension of the sandpile model, 
and certain 3-site probabilities involving a height 1 or a supercritical height. 

Since the bulk of the calculations reported here was done, an article by Jeng has appeared, where 
precisely the same problem is addressed. The two works were carried out independently and differ in 
two ways. The first one is that we study the dissipative model, while Jeng considers it at criticality. 
This fact enables us to deduce the field identifications for the boundary height variables h > 1 from 
the spanning tree characterizations of a single insertion of such variables, technically much simpler to 
solve than the 2-site insertions. The latter are only used as cross-checks, in contrast to the approach 
at criticality which needs them as inputs. So working off criticality offers a simpler and more reliable 
access to the fields. Secondly, we use a different characterization of the height variables bigger or equal 
to 2 in terms of spanning trees, which is based on a modified, two-step burning algorithm. This, we 
believe, leads to a more transparent formalism which generalizes to multisite probabilities. Our results 
and conclusions however fully agree with those of Jeng. 

2 The dissipative sandpile model 

Consider a finite portion L of a square lattice and define at each site i a (sand) height variable hi which 
can take the integer values 1,2,3,.... A configuration C of the sandpile is the set of height values {hi} 
for all sites. The dynamics is defined in terms of a symmetric toppling matrix A. Its entries are all 
integers, positive on the diagonal, negative off the diagonal, and it has row sums which are non negative. 
A configuration is called stable if all heights satisfy hi < An. 

The system evolves in discrete time as follows. To the stable configuration Ct at time t we add a sand 
grain at a random site i (chosen with uniform distribution say), namely we set hi — > /i^ + 1. This new 
configuration, if stable, defines Ct+i- If it is not stable, the unstable site i topples: it loses A^^ grains. 



2 



every other site j receives — grains, whereas X^jgL ^^^^ grains are dissipated (they fall off the 
pile, to a sink). That is, when a site i topples, we update the heights according to 

hj — > hj — Aij , Vj e L. (1) 

If other sites become unstable after the toppling of the site i, they topple following the same rule. All 
unstable sites are then toppled until the configuration becomes stable again. This configuration is then 
taken as Ct+i- In this way, the toppling at the seeded site can trigger a potentially large avalanche, 
resulting in a configuration Ct+i which can be completely different from Ct- 

Provided there are dissipative sites, i.e. sites k for which J2jeL ^'^i ^ ^' ^^"^ dynamics is well defined: 
it does not depend on the order in which the sites are toppled (the model is Abelian) , and the relaxation 
of the seeded configuration to Ct+i requires a finite number of topplings. 

One is generally interested in the long time behaviour of the sandpile. As shown by Dhar [H], 
this behaviour is characterized by a unique time invariant probability measure P£, which specifies 
the probabilities of occurrence of all stable configurations, independently of the initial configuration. 
The moments of this measure, in the thermodynamic limit |L| — s- oo, are what we want to put in 
correspondence with the correlators of a confornial field theory. 

When the dynamics is started from a certain initial configuration, it produces at later times two 
kinds of configurations, called transient and recurrent in the terminology of Markov processes. The 
transient configurations are those which occur a finite number of times only (they may not occur at all, 
depending on the initial configuration). In the long run, they are not in the image of the dynamics, and 
have a zero measure with respect to P£. A simple example is the configuration with all hi = 1, but 
more generally, any configuration with two I's at neighbour sites (sites j with Ay ^ 0) is transient. 

The non transient configurations are recurrent. Their number is equal to det A, the determinant of 
the toppling matrix, and asymptotically, they occur with equal probability, so that the measure P£ is 
uniform on them 8 . A criterion to decide whether a given configuration is recurrent or transient is 
based on the notion of forbidden sub-configuration (FSC): a sub-configuration, with support /C C L, is 
said to be forbidden ii hi < — X]j"eA:\{i} ^3i^ sites i of /C. For instance two neighbour sites with 

heights 1 form an FSC. Then a configuration is recurrent if and only if it contains no FSC 

A practical way to test a configuration is to use the burning algorithm jHj. At time 0, all sites in 
L are unburnt and we define /Co = L to be the set of unburnt sites at time 0. The sites i of /Cq such 
that hi > — X]jeK;o\{i} burnable at time 1. So we burn them, obtaining a smaller set /Ci C /Cq 

of unburnt sites at time 1. Then the sites of /Ci which are burnable at time 2, i.e. those satisfying 
hi > — X]j"gA:i\{i} ^jii burnt. This leaves a smaller set IC2 C /Ci of unburnt sites at time 2. This 
burning process is carried on until no more site is burnable, which means that /Ct+i = /Ct for a certain 
T. Then the configuration is recurrent if and only if all sites of L have been burnt (/Ct = 0)- Otherwise 
/Ct is an FSC. 

The burning algorithm allows to define a unique rooted spanning tree on a graph L*, from the path 
followed by the fire in the lattice [H]. The graph L* has the sites of L and the sink as vertices, and 
has links defined by A: an off-diagonal entry A^- = —n means there are n bonds connecting the sites 
i and j, and each site i is connected to the sink by a number of bonds equal to X^jeL — 
number of grains dissipated when i topples. At time 0, the sink is the only burnt site and forms the 
root of the tree. In the next steps, the fire propagates from the sink to those sites which are burnable 
at time 1, then from the sites which have been burnt at time 1 to those which are burnable at time 2, 
and so on. If a site burns at time t, one says that it catches fire from one among its neighbours that 
were burnt at time t — 1 (or from the sink site at time 1). In case there are more than one of these, a 
fixed ordering prescription is used to decide along which bond the fire actually propagates. (The precise 
prescription will not be needed in what follows; the interested reader is refered to [H] for an example of 
such a prescription.) The collection of all bonds forming the fire path defines a spanning tree, rooted in 
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the sink, and growing towards the interior of the lattice L. 

This improved algorithm establishes a correspondence between the set of recurrent configurations 
on L and the set of rooted spanning trees on L*. The precise mapping, although one-to-one, however 
depends on the prescription used. The prescription we will use below differs slightly from the one defined 
in j^, but is equally valid. The specific sandpile model we consider in the next sections is defined on 
the discrete upper half plane L = Z> x Z, and has the massive discrete Laplacian as toppling matrix 
subjected to the two different boundary conditions "open" and "closed" on the boundary, which we take 
to be the line y = 1. The two toppling matrices are almost identical and differ only along the boundary. 
They both depend on a positive parameter i, controlling the rate at which sand is dissipated when sites 
topple. They read explicitely 




ii i — i are off boundary, 
ii i — i are on boundary. 



A°P = <; -1 if z and J are n.n., A^j = <( " ^ ^ ' (2) 

' ;i • — 1 if z and 7 are n.n., 

otherwise, | 

otherwise. 

One easily sees that t grains of sand are dissipated (transfered to the sink) each time a site topples, or 
< + 1 if it is an open boundary site that topples. This model will be called the massive Abelian sandpile 
model (MASM) in reference to the massive discrete laplacian where ^/i plays the role of a mass. The 
usual, critical model originally defined in j^, is recovered at t = 0. 

In terms of the graph L*, bulk sites and closed boundary sites have a t-fold connection to the sink, 
while open boundary sites have a (i + l)-fold connection to it. In addition, all nearest neighbour (n.n.) 
sites on L are connected by a single bond. 

An easy corollary of the above burning algorithm is that a site with height smaller or equal to the 
number of its neighbours on L (3 or 4) is never burnt at time 1, and therefore catches fire from one 
of its neighbours and not from the sink. Conversely, a site with height strictly larger than A^^ {t ~ 0) 
— which we call a supercritical height — , or an open boundary site with ft, = 4 is set afire by the sink. 
Supercritical height values are those which exist only when t is non zero. 

According to the definition of A, t should take integer values. However the MASM correlations decay 
exponentially, with a correlation length that diverges only when t goes to 0, like The large 

distance limit of the lattice correlations must therefore be accompanied by a small t limit, in such a way 
that their scaling limit ^/t =^ Ma Q,\i — j\ = \z\/a oo he well defined when a ^ 0. So in practice, 
one expands the lattice MASM correlations in powers of i, and selects the dominant terms. These define 
correlators of a massive field theory, which, in this case, turns out to be a massive perturbation "B* of 
c = — 2 logarithmic conformal theory [111 1121 [T5| . 

J {doBe + Ai^ee/A), (3) 

where 9, are anticommuting scalar fields. 

In the course of the calculations, we will make an extensive use of the inverse toppling matrix A~^. 
As is well-known, the inverse of A on the upper half plane can be obtained in terms of the inverse 
massive Laplacian on the full plane A~^, via the image method. For i = (mi, ni) and j = (m,2, 712), the 
explicit formulae read 

(A°P)^i = A^i - A^i = A^i - A-i , f = (m2, -n2), 

(A^i)-.i = A^r.i + A^r.i = A.r.i + A-v^. , ^ (m^, 1 - n^). (4) 

The horizontal translation invariance is preserved in both cases, so that the entries of the inverse matrices 
depend on |mi — TO2I, «i and 712- A short review on values of the inverse massive Laplacian on the plane 
can be found in 

The lattice open boundary condition is identified with the Dirichlet condition in the continuum 
(0 = = on R), whereas the closed boundary condition corresponds to the Neumann condition 
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{do — 89 = dd — 86 — on M.) . The Lagrangian ^ then imphes the foUowing Green functions on the 
upper half plane 

{e{z)e{w)) - {9{z)9{w)) = 0, (5) 
{e{z)9{w)),p = Ko{M\z - w\) - Ko{M\z - w\), (6) 
{e{z)9{w)),, = Ko{M\z - w\) + Ko{M\z ~ w\), (7) 

where Kq is the modified Bessel function. 

3 Unit height and supercritical height variables 

Multisite probabilities for a number of sites to have height equal to 1 or supercritical height values 
{hi > An — t) is fairly easy if one uses the Bombay trick, a beautiful technique designed by Majumdar 
and Dhar 3 . It can be formulated in terms of height configurations or in terms of spanning trees. In 
this section, we will use it in terms of heights, the formulation with trees being a particular case of the 
general characterization given in the next section. 

Suppose that we first want to compute the probability P[hi„ = 1] that a certain site iq has a height 
equal to 1. That probability is simply equal to the number of recurrent configurations with a height 1 
at io divided by the total number of recurrent configurations, which we know equals det A. 

The idea of is to define a new sandpile model in which the height at io is always 1, and such that 
any recurrent configuration of this new model is in correspondence with a recurrent configuration of the 
original model where the height at iq is 1. To freeze the height at io to the value 1, one simply reduces 
the diagonal entry of the toppling matrix to 1. So the toppling matrix A' of the new model will have 
^'ioio ^ ^' Then in the new model, the site zq will topple whenever its height exceeds 1, and each time 
it topples, it will lose a single grain which will go to one the neighbours or to the sink. Consequently, 
?o will have a single connection, either to the sink or to one its neighbours in L. Finally, the neighbours 
of io cannot have a height equal to 1 in a recurrent configuration, so that they assume only A^^ — 1 
values. This can also be enforced in the new model by decreasing the diagonal entries of A by 1 for 
those neighbours of io which are no longer connected to io- As the connections fix the off-diagonal part 
of the toppling matrix, and the height ranges fix its diagonal part, this will determine A'. 

Thus the number of recurrent configurations with a height 1 at io is equal to the total number of 
recurrent configurations of the new model, itself equal to the determinant of the new toppling matrix 
A'. Setting A' = A + B'-'o\ one obtains 

^^'"'^ = ^] ^ ^ '^"^^^ + A-ii?(^°)), (8) 

where A is the toppling matrix appropriate to the boundary condition one considers. Because the 
difference A' — A = _b(*o) is non zero only on sites around io, the previous formula reduces to the 
calculation of a finite determinant, even on an unbounded lattice L. 

On the discrete upper half plane, the defect matrix B^'o) depends on the location of io- If io is off 
the boundary, and if one keeps it connected to one of its four neighbours, then _B(*") is equal to 
/-3-t 1 1 1 \ 

10-10 ^ ^ 

V 1 0-1/ 

on io (first label) and any three neighbours of io, and is identically zero elsewhere (if the only connection 
of io is to the sink, the B matrix is 5-by-5). In this case, the probability is given by a 4-by-4 determinant, 
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and depends on the distance m of iq to the boundary. At the critical point [t — 0) and for large values 
of TO, it is equal to ^1] 

F[;i,„ = i] = Pi[i±^ + ...], (10) 

where the + (resp. — ) sign refers to the open (resp. closed) boundary condition. Pi = 2(7r — 2)/-k^ ~ 
0.0736 is the probability that a site deep inside the lattice (equivalently, on the infinite plane) has height 
1 in the critical ASM. 

If iQ lies on the boundary, the matrix B^'o) depends on the boundary condition and takes one of the 
two forms 

/-i-t 1 1 \ /-2-t 1 1 \ 

B^o) =1 -10 , S^i"^ =1 -10 . (11) 
V 1 -1/ V 1 -1/ 

The corresponding critical probability P[h^-, = 1] is then a constant, which only depends on the type of 
boundary the site iq is on 

pop _ 9 _ 42 320 _ 512 pel _ 3 _ 2 
1^2 ^ + 3^2 9^3' ^1-4 ^- ^^^) 

The one-site probability P[hif, = 1] can easily be computed for t 7^ 0, but will not be needed in what 

follows. 

The probability that a site be supercritical can be treated in a similar way, and is actually simpler. 
One now takes t ^ 0, since the probability does not make sense at i = 0. 

Any site io, whatever its location and whatever the boundary condition, has t possible supercritical 
height values, namely h = Ai^i^ — t+l, . . . , Ai^ig . The probability that a site io has a fixed supercritical 
height h does not depend on h, because a recurrent configuration remains recurrent if one replaces a 
supercritical height at io by another one. Therefore one has 

P[io is supercritical] — P[hig > Aigig — t] — t P[hi^ ~ h], (13) 

where h is any fixed supercritical value. 

It is actually easier to compute the probability that zq is not supercritical. To do that, one has to 
count the recurrent configurations with hi^ < Ai^ig — t (= 4, or 3 on a closed boundary). In a new 
model defined by the new toppling matrix A'^j — Aij — t Si^ig Sj^ig on the same lattice, all recurrent 
configurations have iq not supercritical. One obtains 

P[io is supercritical] = 1 - = t A.^J^ , (14) 

and 

P[h,g ^ h] ^ A-^. (15) 

The defect matrix method works here too, and is simpler because the appropriate matrix 3^^°^ — 
—tSi^ig 6j^ig has rank 1. The corresponding one-site probability is then given by a 1-by-l determinant. 

Let us note that the probabilities H15|l are well-defined for any strictly positive value of t, but behave 
badly in the critical limit i — > 0. For a closed boundary condition, they have a logarithmic singularity at 
t = 0. For an open boundary condition, they have a finite limit at i = 0, but which is not a probability 
in general: A^^^ — \ — 2/it = 0.3634 for iq on the boundary, and then grows logarithmically with the 
distance of «o to the boundary. 

So instead we will consider the probability H14|l for a site or a collection of sites to be supercritical 
without specifying the actual heights. As we will see below, that observable has well-defined correlations 
in the massive scaling limit, and corresponds to a field that vanishes in the critical limit. 

After the one-site probabilities, multisite probabilities and correlations can be computed by the 
same method almost routinely. The observables we consider in this section are the two boundary 
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°^^'''\x) = lim [S{K/a - 1) - Pi] , (t)'^'''\x) = lim [S{K/a is supercritical) - A A], (16) 



random variables 6(hi — 1) and d{i is supercritical) corresponding to the events "i has height 1" or "i is 
supercritical" for i a site on the boundary of the upper half plane. In order to get fields whose expectation 
value vanishes infinitely far from the boundary, one considers the random variables subtracted by their 
average value. Anticipating the scaling dimension 2 or 4 of the above random variables, we define their 
scaling fields by 

— n '~2 I'^i^^/a - 1) - Pi] > = lim 

subjected to the scaling relations t = a'^M'^ and = Mx. 

To compute n-site probabilities, one simply inserts the proper defect matrices at the locations of the 
observables, so that the full defect matrix is a direct sum of n matrices B^'-* or S'^*^. One should however 
remember that S^^^ is not the defect matrix for i being supercritical but for i not being supercritical, 
the complementary event. The scaling limit of the latter gives rise to a field (/><, from which 0> = —4>< 
is recovered. The probability then reduces to the calculation of a finite determinant whose entries are 
combinations of entries of the inverse toppling matrix. As the scaling limit takes t to zero, one expands 
these entries in power series of keeping only the dominant term. ^From (|16|l the latter yields the field 
theoretic correlation of fields (pi and (/>> , which are then identified with explicit fields of the Lagrangian 
theory Q. This is a main advantage of working with the massive theory that this identification is 
essentially unambiguous. 

The simplest way to proceed to the identification of the boundary fields <f>i and </>> is to use other 
lattice observables with already known field identifications. Examples of such observables are precisely 
the bulk version of the above two random variables. The corresponding bulk fields have been identified 
in (see also ^Hj for a proof that these identifications are consistent with a broad class of multisite 
correlations) , 

_~ _ - A'f2 ^ 

(j)i{z) = -Pi[:d0dd + ded9:+— -.ee-.], (17) 

<P>iz) ^ ^ -.ee-.. (18) 

One first computes the 2-point correlations involving one boundary observable and one bulk observable. 
From them, one may infer what the boundary fields must be, and then cross-check their form from other 
correlations. 

We do not give much detail as the calculations are fairly straighforward, but simply illustrate the 
method in a simple case, namely the identification of 4>'^{x) on an open boundary. 

We take two reference sites, io on the boundary and jo in the bulk, far from the boundary. The 
probability that they both be non supercritical reduces to a rank 2 determinant 

det[A°P + SM + 5O0)] _^ Jl-t (A°P)-i^ ~t (A°P),r^ 



PHoJo non supercrit.l = = det ,\ 'Y" \ . • (19) 

Subtracting the product of disconnected probabilities obtained from (|14|l . one has 

P[ioJo non supercritical] connected = -t'^ (A°p)^:^^-^. (20) 

If one chooses the two sites on a vertical line io — (0,1) and jo — (0,m.), then (see for instance 
Appendix A of 6 ) 

i^n-)o - ^io.imm) - ^ioUmm) - ^^Kim^t) + ■ ■■ (21) 

where the dots stand for subdominant terms in t. The dominant term in the connected 2-site probability 
is thus 

^3 

P[io, jo non supercritical]conncctcd = n K'^{m\/t) + . . . (22) 
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Using the scaling relation H16() . valid also for bulk fields (with a power a^^ for (f)i and (/><, 0>), one 
finds the boundary/bulk 2-point function 

/\,f6 

(0°<P(x) 0<(x + ly)) = - — K',\My). (23) 

^From the explicit form of (l)<{z) — — '-09: given above, one eventually arrives at 

o /v f 4 

(0°P(a;) -M: {x + ly)) = K',^{My). (24) 

— TT 

Using the Green function (jSJ on the upper half plane with an open boundary, one sees that the only 
possible field assignment is <ti°y{x) = 2^ ■.8989:, and therefore ip'^ix) = — -.8989:. 

Proceeding in the same way for the other observables for the two boundary conditions, we find the 
following scaling fields 

/6 _ 150 ^ mA^ ^ 

\7r 37r^ Ott-^ / tt 

<^f ^ :d9d9: +^M^ -M: ] , = ^ -M: . (26) 

TT V4 TT/ '■ 16 27r 

These field identifications have been checked to be consistent with many multisite probabilities: 

2-site and 3-site boundary/boundary probabilities as well as 2-site and 3-site mixed boundary/bulk 

probabilities. 

The massless limit is simply given by the limit t ^ in the MASM and M ^ in the field theory. In 
this limit, the fields 4>> and (jf^ are obviously null and the unit height fields for the two boundaries are 
identical up to a numerical factor. One may note that the latters are proportional to the holomorphic 
stress-energy tensor, and being descendants of the identity, they belong to a chiral representation Vq 
|13j . This is consistent with the fact that the only fields living on an open boundary are fields of 
Vo, and that those living on a closed boundary belong to an TZ^ representation, which contains Vq as 
subrepresentation [TOIITt] . 



4 Spanning tree representation of recurrent configurations 

For the other height variables, the situation is not as easy. Although having a height 1 or a height 
2 at a given site does not seem to make much difference, the counting of the corresponding recurrent 
configurations is technically much more complicated for a height 2 (or 3 or 4) than for a height 1. The 
defect matrix method no longer works^, and the only practical alternative seems to be the use of spanning 
trees. One then clearly sees the difference: in terms of spanning trees, a height 1 is characterized by 
a local property of the tree around the reference site, while the other heights are characterized by non 
local properties of trees 0j. 

As mentioned earlier, the rooted spanning tree is defined on L*, the lattice L augmented by the sink 
site, at which the tree is rooted. All sites i of L are connected by A^i bonds to other sites: —Ay bond(s) 
connecting i to j and X^jeL bond(s) connecting i to the root. With these definitions, the Kirchhoff 
theorem asserts that the number of rooted spanning trees on the graph L* defined by the matrix A is 
equal io N — det A, precisely the number of recurrent configurations. 

As a rooted spanning tree is a connected graph containing no loop, every site i is connected to the 
root by a unique path. A site j is said to be a predecessor of the site i if the path from j to the root 
passes through i, or equivalently, if j lies on a branch growing from i. A site i which has no predecessor 
is called a leaf (the end of a branch) . 

Priezzhev first and then Ivashkevich used the correspondence between recurrent configurations and 
spanning trees to compute respectively the 1-site probabilities in the plane, and the 1- and 2-site 

^Except a lieiglit 4 on an open boundary, which can be handled like a supercritical height. 
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probabilities on the boundary, open or closed, of the upper half plane. For the 1-site probabilities at 
io, they decomposed the set of recurrent configurations into subsets Sa, where Sa contain the recurrent 
configurations which remain allowed for any heights hi^ > a and which are forbidden otherwise. These 
subsets Sa can be characterized in terms of rooted spanning trees and their cardinal can be computed 
using classical results in graph theory. As we will see, they decomposed local tree diagrams as sums of 
non local diagrams. This system is invertible without further input for on a boundary, but is not for 
?o in the bulk. So the calculation of probabilities for sites in the bulk is more complicated. 

For the 2-site probabilities, Ivashkevich used a similar decomposition of the recurrent configurations 
into subsets Sab- This decomposition however raises certain questions, and will not be used here. Instead 
we set up a particular one-to-one map between recurrent configurations and the rooted spanning trees, 
based on the burning algorithm. In the case of single site probabilities, the mapping yields the same 
characterization in terms of trees as the Sa decomposition but is much more transparent in the case of 
multisite probabilities. 

The burning algorithm, as we described it Section 2 complemented with an ordering prescription, 
establishes a one-to-one mapping but with no clear correspondence between the height values at the 
reference sites and the bond arrangements of the trees around those sites (except for supercritical heights 
which are directly connected to the root of the tree). For example, depending on the recurrent configu- 
ration, a site with a height 4 can be a leaf on the tree or can support 1, 2 or 3 branches. To avoid this 
problem, we proceed in two steps as follows, assuming that none of the reference sites is supercritical. 

First, we run the burning algorithm and let the fire propagate through the lattice until no more 
site is burnable but preventing the reference sites from burning. When this is done, one is left with a 
sublattice Lb of burnt sites and a complementary sublattice Lu of unburnt sites. The algorithm, using 
for example the ordering prescription of [HI, yields at this stage the part of the spanning tree on Lb. The 
other part Lu is eventually burnable too and is actually burnt in the second step. By definition, none of 
the reference sites is burnt yet, and a certain number of them, at least one, are burnable. Those which 
are burnable are burnt simultaneously, and trigger the fire propagation through Lu, thereby completing 
the spanning tree to the whole lattice. So the complete tree is made of two pieces, a subtree 7b on Lb, 
and another 7^ on Lu. The subtree 7^ itself may have several roots which are among the reference sites 
those which were burnable and which set fire to the whole of Lu. It is at those sites that is grafted 
to 7b to make the full tree T. As we will see, only the shape of the unburnt sublattice Lu is used to 
characterize the height of the reference sites. 

This slightly modified burning algorithm establishes a well defined correspondence between spanning 
trees and heights of the reference sites in the critical as well as in the massive Abelian sandpile. Let 
us see how this works for the single-site probabilities, and how it allows to compute the 2-site height 
correlations where one the two heights is equal to 1 or is supercritical. At this stage we will be able to 
identify the boundary fields corresponding to all heights. In the following section, we will compute other 
2-site and 3-site correlations to confirm these identifications. 

Let us consider a configuration of the MASM on a square lattice L, and let us focus on a fixed site 
io, the reference site. We will take L to be the upper half plane, but what follows applies to any sort of 
portion of Z^, bounded or unbomided. 

• If the height at is supercritical, then is set afire by the root (is burnt at time 1). Thus all trees 
corresponding to those configurations have a bond connecting the root and ig. The probability 
that io be supercritical is thus 

P[io is supercritical] = (27) 

where JV — det A and AC.if, is the number of different spanning trees which "use" one of the t 
bonds between the root and ig. One way to compute A/'^^i^ is to modify the toppling matrix by 
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removing the connections between io and its nearest neighbors on L so that iq has connections only 
to the root. Then Af^So = det A' with A' = A + S''^*"' and the finite-dimensional defect matrix 
given by S-^-] —t — ^igig, 5',^°] = sf°^^ ~ I for the nearest neighbors of zq, and zero elsewhere. 
A simpler way is however to compute M — the number of trees which do not use the bonds 

between io and the root. This can be done by removing precisely these bonds and leads to the 
1-dimensional defect matrix 3'--'^°^ of the previous section, with the result given in H14|) . 
The same arguments apply to a height 4 at an open boundary site. The only difference is that 
such sites have t + 1 connections to the root, so that 

^4 = J ^^^^p- = ' *o ^P^"^ boundary. (28) 

If the height at ig is less or equal to the number of nearest neighbours on L, we use the burning 
algorithm to define a partition L = Lb U Lu as explained above. As one looks here at 1-site 
probabilities, there is only one reference site, so that in the sublattice — L^q , io is the only 
burnable site after all sites of Lb have been burnt. It is therefore the root of the subtree 7u = 7^^ . 
The height at iq can now be related to the properties of the subtree . 

If Li„ contains only the site iq, then hig can take any of the values 1,2,..., rii^ , where n^^ = 3 or 4 
is the number of nearest neighbours of iq. The full tree T is simply obtained by connecting io to 
7b through one of the rii^ bonds connecting iq to its nearest neighbours, so that io is a leaf on T 
(io must be connected to a nearest neighbour in L and not to the root, since it catches fire from 
one of them) . If one denotes by Afi the number of spanning trees on L where io is a leaf grown on 
one its neighbours, then 

P. = (29) 

If Lig contains one nearest neighbor of the site io, the value of the height at io must be in the set 
{2, . . .,nig} (it cannot be 1 since otherwise io would not be burnable and could not set L^^ afire). 
There are now riig — 1 possibilities to connect 7^^ to Tb, one for each nearest neighbour of io in Lb- 
They correspond to the height values compatible with the burning algorithm and thus the height 
2 probability reads 

P..^ + ^_^. (30, 



N (n,o -l)A^' 



where M2 is the number of spanning trees where io has exactly one predecessor among its nearest 
neighbours. 

The higher height probabilities can be determined by the same arguments, 

^- = ^--^+ (.,-t+l)AA ' ^ = 1'--'-- (31) 

where A/fc is the number of spanning trees T on L in which io has exactly fc — 1 predecessors among 
his neighbours. 

One sees that the computation of the various 1-site probabilities requires the calculation of the 
numbers Mk- If the lattice L is the discrete upper half plane, and for io a site on the boundary, the 
Figure 1 describes the types of trees which contribute to the different NkS. 

These diagrams represent the restriction of trees to four sites, namely the reference site io and its 
three neighbours (the dashed line represents the border of L, pictured as the lower half plane !). The 
arrows indicate the direction of the path towards the root (opposite to the fire propagation line). The 
black dots are the nearest neighbours which are predecessors of io, the white dots are those which are not. 
Those diagrams labelled by identical greek letters contribute the same amount to the corresponding A4- 
The tilded letters refer to diagrams which are the mirror images of the diagrams with untilded letters. 
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M = La + L±J + 

ai a2 0(3 



A/-2 = tr:] + m + o + u:] + Li:] + m 

Pi P2 h h 71 72 




Figure 1: Non local diagrams contributing to the 1-site probabilities. 

The diagrams in Figure 1 represent non local constraints on the compatible trees. The presence of an 
arrow between io ^md one or more neighbours poses no computational problem, because it only means 
that the tree has to use specific bonds, and the counting of those trees can be achieved by modifying 
locally the toppling matrix by a finite rank defect matrix. But for a neighbour to be a predecessor or 
not is a non local property, and enumerating the relevant trees is trickier. For the trees contributing to 
As for instance, one sees that a nearest neighbour, call it ii, of zq can be a predecessor of zq because ii 
is connected to io through the nearest neighbour bond {ii catches fire directly from zq, like in (5, e and 
e), or through a long path around the lattice (zi catches fire after a long sequence of burnings, long after 
zo got burnt, like in (j) and 0). 

Various classes of non local diagrams in Figure 1, with their constraints on predecessorships, sum up 
to local diagrams, where these constraint are relaxed and only local arrow configurations are imposed. 
As the local diagrams are easily calculable by toppling matrix adjustments, this yields linear relations 
among the non local contributions. For a generic position of the reference site iq, the linear system is 
not invertible and is not sufficient to compute the non local contributions. The crucial observation made 
by Ivashkevich in [S] was to note that it actually becomes invertible if zq is on a boundary, which allows 
to reduce the non local contributions to local calculations (this statement will have to be qualified in the 
case of 2-site insertions). Indeed the relations read explicitely |5j (as we consider the heights 1, 2 and 
3 only, we do not need to distinguish the diagrams for closed and open boundaries, which only make a 
difference through the proper toppling matrix to be used in the explicit local computations; at the same 
time, that also makes the number of equations smaller) 

Qo - lZ] = Cll 

a 





e e 

The four equations on the right are clearly obtained from the corresponding four on the left by a 
mirror symmetry about the reference site zq. As the toppling matrix is invariant under that symmetry. 
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the three tilded diagrams contribute the same amount as the untilded ones, f3 = (3,e = e,(l) — (j), as do 
the tilded and the untilded local diagrams. Thus the four equations on the right are redundant, and 
one is left with the linear system on the left. This system is manifestly invertible for the six non local 
contributions, noted a, (3, 7, 5, e and </>. 

Let us note that the non local diagrams a, (3, 6 and e turn out to be entirely local, because the 
arrow configurations make the predecessorship properties redundant. As the height 1 probability Pi is 
given solely in terms of a, its computation is purely local. This remains true for any multisite height 1 
probabilities and for arbitrary positions, in the bulk or on boundaries. 

The Kirchoff theorem allows the local diagrams to be computed by the defect matrix method. The 
presence resp. the absence of an arrow from i to j means that one counts all trees which contain resp. 
do not contain that oriented bond. One defines a new toppling matrix A' by setting to — e the i,j entry 
if the i j bond is to be used in the tree, and to if that bond is not to be used; moreover the diagonal 
entries of A'^^ = An should remain equal to the number of bonds going out from i. Then if n bonds are 
to be used, the determinant of A' will contain a highest degree term e" whose coefficient is the number 
of trees which precisely use the given n bonds in the prescribed direction (see for instance 0]). 

Writing as before A' = A + B, one finds for Q4 for example 

/-3+e 1 1-e 1 \ 



B 



1-e -1 + e 
1 0-10 
V 1-e -1 + eJ 



(33) 



and, in the critical limit, 

n. A 1 f i — 4 on closed boundary, 

^ = 77 = lini 4 det[I + A-'B] - \ " 3 (34) 
J\l J\l e->oo e ^ { -^^^^ ) Qj-^ open boundary. 

The calculation of the other five local diagrams and then the inversion of the linear system yields the 
values of A/i,A2 and As, and in turn of Pi,P2 and P3. In the critical limit, one recovers the numbers 
given in |S]. 

In order to identify the height boundary fields, we need 2~site correlations involving the boundary 
heights 2 and 3. Again the simplest is to look at the correlations of a boundary height 2 or 3 with 
a known boundary variable, namely a height 1 or a supercritical height value. The advantage is that 
the latters are already known from the previous section, but more importantly, they correspond to 
local defect matrix insertions. This makes the above formalism, useful to compute 1-site probabilities, 
essentially valid. 

Because one can force a site iq to have height 1 or to be supercritical by modifying the toppling 
matrix by A — > A(io) = A + or A + S"''") (or better A + S'^")), the 2-site probabilities P[hi„ = 
1 or supercr., = 2 or 3] can be viewed as 1-site probabilities for the height at but with the 
toppling matrix A(«o) to account for the constraint at Iq. Then the above method remains completely 
valid provided we replace A by the appropriate A(io), itself to be modified by matrices like in (|33|l in 
order to compute local diagrams ^ . If one does that and use A(io) as the normalizing toppling matrix, 
one is actually computing the conditional probability for having a 2 or a 3 at jo conditioned on having 
a height 1 or a supercritical height at iq. To get the joint probabilities, one simply multiplies the final 
answers by P[hi„ — 1] or P[hi„ — supercr.]. 

The non local diagrams contributing to the numbers A/l- remain as in Figure 1. However the tilded 
and the untilded diagrams no longer contribute equally because the mirror image about jo spoils the 
constraint at iq, and does not leave the toppling matrix A(io) invariant. Therefore it is the full system 



^To keep the decompositions of the Mk in terms of the non local diagrams as in Figure 1, the site io should not be too close 
to jo. 
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(|32|l that needs be solved. It is overdetermined as it involves 10 equations for only 9 unknowns, but the 
number of equations is reduced by one due to the following identity 

P(l) + P(2) + P(3) = P(l) + P{2) + P(3), (35) 

satisfied for all values of i as a simple consequence of the fact that the inverse of A(io) satisfies a discrete 
Poisson equation. The procedure is otherwise identical to that for the 1-site probabilities. 

For the open boundary, the boundary joint probabilities of a site with height 2, 3 or 4 and a site 
with height 1 or with a supercritical height all have the same form as two unit heights on the boundary. 
It means the same field identification up to a numerical factor: 

= {-'4 + W^-f^)--ded~9:, (36) 
'^T - C4-^ + m--dOd9:, (37) 



bT - 



= -^-.dedO:. (38) 



The last line can most easily obtained from (|28() . which implies 

For the closed boundary, the correlations involving a height 2 or 3 have a more complicated structure. 
For example, one finds (with m = \iQ — jqI) 

P[h, > 3, hj, - 2] - 1^ [^K^imVi) - (3 - ^)K'i{m^t) - ^K^{7n^t)K'^{7n^t)\ + ... (39) 
P[h, > 3, hj, = 3] = 1^ [\Ki{mVi) - ^K'^\mVt) + \KQ{'m^ft)K'^{mVt)\ + ... (40) 

These results and the corresponding ones for hi^ — 1 are compatible with the following field assignments 
for the height 2 and 3 



' = (f - i) --dOd-e: +^ -.eddO: +{^ - ^)KP -M:, (41) 
- 8 -dOdO:-^ -.edde-.-^KP -.eO: . (42) 



V^2 

This identification is in this case not unique, since the field theory is invariant under 9—^9 and 9 —^ —9. 
One could in particular change : 9{dd9) : for : {dd9)9 :, which are different fields since their correlation 
contains a logarithm while their self-correlations do not. If one requires that the sum (pf + (f)f + (pf be 
zero in the massless limit, then the choice for one the two fields must be made for both (pf and (pf. The 
correlations computed in the next section confirm this. Note that the sum of fields (pf + (p2 + (p'i + (py 
vanishes identically. The similar sum in the open case, cp'i' + (p'2' + + ^4'' + (p'y' , vanishes at the 
critical point only, because the dimension of does not match the dimension of the universal terms 
of the other fields. 

At this stage, all boundary height fields for the massive Abelian sandpile model have been determined. 
To have more checks on the field identifications, we compute in the following section, all 2-site and some 
3-site height correlations. 



5 Higher boundary correlations 

In the previous section we have seen that the multisite probabilities where only one reference site iq 
has a height value in {2, 3} can be computed from the diagrams listed in Figure 1 by using a toppling 
matrix properly decorated by defect matrices to account for height constraints (height 1 or supercritical) 
at the other sites. The calculation of multisite probabilities where two reference sites zq and jo have a 
height value in {2, 3} leads naturally to pairs of such diagrams, one at iq, the other at jo- However the 
situation becomes technically more complex because sites in the diagram at zo can be predecessors of jo 
and/or the other way round. So the topology of the spanning trees can be more complicated and their 
counting more difficult. 
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Let us first consider the 2-site probabilities Pat = P[hig = a, hj^ = b] for a, b in {2, 3} and where 
io and jo are on the boundary of the upper half plane. We start the burning algorithm as explained in 
Section 4 without ever burning the sites ia and jo, and until no other sites than those two are burnable. 
This yields a sublattice Lu = hi^jg of unburnt sites, which subsequently catches fire either from io or 
from jo, or from both if they are both burnable. In turn the fire propagation on Lj^j^ defines a subtree 
7^ = T^ojo , rooted at io, or at jo, or at both sites. The full tree T is made up of the subtree T^, living on 
the sublattice of burnt sites, to which is grafted at io and/or jo- 

The restriction of any tree to the neighbourhood of a reference site looks like one of the non local 
diagrams shown in Figure 1. So one can visualize the restriction to the two neighbourhoods by a pair of 
such diagrams. Using the same labelling as in Figure 1, we will denote the pairs of diagrams by pairs of 
greek letters (with indices), the first one for the diagram around io, the other for the diagram at jo. In 
an obvious notation, a pair of greek letters belongs to a certain set Afk x J^i- As we did in Section 4 for 
the one site probabilities, we have to compute which probabilities Pat a pair of diagrams contributes to. 

For 1-site probabilities, we know from Section 4 that the diagrams in A/fc contribute equally to the 
probabilities Pa for fc < a < 3. Indeed the three diagrams ai, o.^, as of M\ are obtained from each other 
by changing the arrow around the reference site. The change converts a tree which is compatible with 
a diagram on into a tree which is compatible with another diagram a^, and this fact shows that the 
number of trees compatible with a diagram q,; does not depend on i, namely a\ — oii — olz or M\ — 3a. 
As the position of the arrow determines univoquely the height value, the three probabilities P\, Pn, P3 
get an equal contribution MxjZM from the diagrams in M\. The same is true of the six diagrams in Mi. 
They come in pairs {Pi, (Si), {^1,^2), (71,72), where the diagrams within a pair are related by changing 
the direction of the arrow coming out from the reference site. The same arguments as above show that 
/3i = /32, = /32, 71 = 72, and that P2, P3 receive an identical contribution /3 + /? + 7 = M2/2M from 
the diagrams in Mi- For A/3, each diagram is on its own and contributes to P3. 

In the case of 2-site probabilities, the same arguments would show that the diagrams in A4 x Mi 
contribute equally to the probabilities Pab forfc<a<3,/<6<3, provided one can prove that 
changing the direction of an arrow in the way recalled above in either diagram, or in both diagrams, 
turns a compatible tree into a compatible tree of the same class. Because the two diagrams can now be 
linked by fire paths, this is no longer guaranteed, and actually fails in a few cases, pictured in Figure 2. 

On the first line of Figure 2, one sees for instance that the diagram denoted by is a pair 72 /3i. It 
is linked in such a way that when one changes the arrow in j3\, one obtains a well-defined tree (noted 
Ai) compatible with a pair cfibPi- If one changes in the arrow of 72, one obtains the diagram Ai, of 
the type ^i4>a- Changing the arrow of 72 and of j3i introduces a loop, and so cannot contribute to a 
2-site probability. The trees compatible with the diagrams Ax , Ai and A^ are related by local changes 
of arrow, but belong to different classes, namely A/2 x A/3, A/3 x Mi and Mi x Mi- There should normally 
be a fourth diagram, in A3 x M-a , but which does not exist as a tree. 

It is not difficult to see that the misbehaviours with respect to arrow changes can only be of the type 
shown by the triplet {Ax,Ai,Az)- When the two diagrams are tied in a special way by the fire paths, 
one change of arrow in a diagram in Mi x Mi sends it to a diagram in Mi x A/3 or A/3 x A/2, and two 
arrow changes introduce a loop. 

Figure 2 shows four triplets of diagrams where this peculiar behaviour occurs. Diagrams labelled by 
the same capital letter are in equal number, since the numbers of compatible trees are equal. The twelve 
diagrams shown in Figure 2 and the mirror diagrams (not shown in Figure 2), obtained by exchanging 
the diagram at io with the reflected one at jo and vice-versa, make the complete list of misbehaved 
diagrams. We will denote the mirror diagrams with tildes. 

The two-step burning algorithm allows to determine which probability each diagram contributes to. 
In the diagram Ax for instance, the subtree catches fire from the eastern neighbour of jo- It is thus 
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A. 



B^ 




Ci 



Co 



Figure 2: Non local diagrams representing spanning trees which have an anomalous behaviour under a local change 
of arrow around «o and/or jo- The mirror diagrams must be added to have the complete list of such diagrams. 



burnable at a time where only one of its neighbours is burnt, and so must have a height 3. The other 
reference site iq is not burnable at the time Ti^ja catches fire despite the fact that its western neighbour 
was burnt, which implies that its height is at most 2. When iq is burnable, it has two burnt neighbours 
and one southern unburnt neighbour, meaning that its height must be 2. Thus Ai contributes to P23- 
One finds similarly that the first column in Figure 2 contribute to P23, the second column to P32, and 
the last column to P33. 

We define [Mk x A/j] to be the set of trees in Mk x Ni which do not have this sort of misbehaviour 
under a local change of arrow. The set Nk x Ni is equal to [Mk x Ni] except in the following three cases, 

AA2 X M = [N2 X N2] + A3 + S3 + C3 + + i3 + -Ba + C3 + Cg, (43) 
AA2 X AA3 = [N2 X AA3] +Ai+Bi+Ci+Ci + A2 + B2 + C2 + C5, (44) 
AA3 X M = [N-i X N2] +A2 + B2 + C2+ C5 +Ai+Bi+Ci+ Ci. (45) 

The trees in [Mk x Ni] contribute equally to the probabilities Pab^ k < a < 3 and ^ < 6 < 3, while 
those compatible with the diagrams of Figure 2 must be handled separately. One obtains 

[M2 X M2] 



P22 


- P12- 


^P21 


-Pll 


P23 


= Pl3- 


^P22 


-P12 


P32 


= P22 - 


^P31 


-P2I 


P33 


— P23 - 


HP32 


— P22 



4Af 
[M2 X N3] 

2JV 
[Ms X JV2] 

2N 
M3 X 7V3 



[Ai 



f 



[Ai 



A2 + B1 



[A3 



+ B2 

S3^ 



+ C2 

^3 + 



^3 + 



S2 + Ci + C2 + C4 + C5 + mirrors] 



-B2 + C2 
-Bi + C'i 
B3 + CS + 



fC'4] 



(46) 
(47) 
(48) 

(49) 



The subtracted term in P33 is due to the fact that the part of P23, P32 related to the misbehaved diagrams 
in Figure 2 (first and second columns) do not contribute to P33. 

The sets [Afc x A/j] will be further partitioned in classes labelled by a pair of diagrams, f.i. [A/a x A2] — 
[PiPi] + [P1P2] + • . .. One will remember that the cardinal of a class does not depend on the numerical 
indices attached to diagrams, so that = |[/3i/32]| = |[/32/3i]| — |[/?2/32]|, and so on. Replacing 

however (j)a by (j)b in a class does not necessarily conserve the cardinal of that class, f.i. \[P(t>a]\ 7^ |[/3</'b]|- 
The 2-site probabilities can be computed if the numbers of trees in these subclasses and of those 
compatible with the non local diagrams of Figure 2 can be calculated. 

As for the 1-site probabilities, we can decompose each pair of local diagrams as a sum of non local 
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ones. We have for example, 

Qi.i = 




Q2,3 = 



m ( 



(51) 



As such, the hnear system one obtains in this way is underdetermined. Let us proceed to the counting 
in the general case, that is, when heights equal to 1 or supercritical heights are inserted at other places 
than zq and jq. In this situation, the full system must include all pairs of local diagrams (for example, 
one would have the equation Qi,i and its mirror image Qi.i)- 

There are 81 equations like USUI) and (|51|l . since every such equation is a pair of local diagrams, chosen 
from the nine diagrams appearing on the last five lines of Eq. 132|) . There are 9 independent variables 
for the classes of [A/'2 x ^2] (pairs of elements in /3, 7} since the indices are irrelevant), 21 variables 
for the classes of [1^2 x Nz] and [A/a x M2], and 49 variables for A3 x A/3. To these one must add the 
four variables A, B, C and C, for the diagrams of Figure 2 (one can show that A = A and B = B). 

In total, one has a linear system of 81 equations for 104 variables. It is actually worse because the 
equations are not all independent, due to some non trivial identities among local diagrams (like in (|35|) ). 
It is however possible to compute the probabilities in terms of a reduced number of variables. 

In the calculation of the 1-site probabilities, the non local diagrams 4>a and (j)h each brought an equal 
contribution, because for any tree compatible with (pa, there is a tree compatible with 0;, and vice-versa. 
Thus a single variable (f) was used for the two diagrams. 

The substitution of (pa by 0t, in a pair of diagrams does not always conserve the number of trees, 
so that the number of independent variables for pairs of diagrams involving a (p cannot be reduced by 
a factor 2. However, one may separate in [M2 x A/s], [Nz x N2] and A/s x A/s the trees for which the 
substitution is allowed from the others, like what we did above regarding the change of arrows. 

It turns out that this is useful because only a reduced number of pairs of diagrams misbehave under 
the change (pa ^ (pb- Up to mirror symmetry, they are all given in Figure 3. 

For instance the two diagrams Ai and A2 are pairs [/920a] and [/3i(/)a], contained in the set [N2 x Nz\- 
The change (pa — *■ (pb in Ai requires a change of direction in the path going from the southern neighbour 
of jo through and back to joj which is not possible. The pair Hi and 112 corresponds to the diagrams 
[(3i(pb] and [(32<Pb], also in [A2 x A3]. Their permuted versions belong to [A3 x 7V2], while all the other 
diagrams in Figure 3 arc in A3 x A/3. The diagrams whose labels differ only by the numerical subscript 
contribute equally, so Ai — A2 but Ai ^ Ai. 

If one denotes by curly brackets the sets of trees which are closed under the change (pa <-> (pb, one 
can write 



W2 X AA3] = {[AA2 X AA3]} + Ai + A2 + Hi + n2 + A'' 4 
[AA3XAA2] = {[AA3 X AA2]} + A? + A^ + nf + n^ + AH 

AA3 X 7V3 = {7V3 X 7V3} + A3 + A4 + n3 + n4 + $1 + 

+ mirrors and permuted 



- A2 + ni + n2, 

$2 + ^'l + *2 + ^1 



(52) 
(53) 

(54) 



The diagrams Ai + A2, Hi + 112 and the permuted tilded versions contribute equally to P23 and 
P33, since they are in [1^2 x A3]. The diagrams A^" + Aj, II^" + 112 (and the permuted tilded versions) 
contribute equally to P32 and P33, whereas all the others contribute to P33 only. Thus the expressions 
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no 





Figure 3: Non local diagrams representing spanning trees which have an anomalous behaviour under the substitution 
of 4>a by or of (pi, by (pa- The superscript p indicates that the two diagrams at the reference sites have been 
permuted. All mirror diagrams must be added. 



for the 2-site probabilities become 

[M2 X A/'s] 



P22 
P23 



P12 + P21 — Pii 
P13 + P22 — P12 



{[AA2 XAA3]} 



P 



33 



2Af 2Af 
Ai + A2 + Hi + n2 + A? + + + ftP] 
X Af,} 



(55) 

[2Ai + 2Bi + 2C\ + 2C4 + 2A2 + 2B2 + 2C2 + 2C5 

(56) 



P23 + P32 
1 



Pr. 



N 
1 

N 



[A3 + B3 + C3 + Ce - - A2 - Si - B2 - Ci - C2 - C4 - C5 + mirrors] 
[A3 + A4 + 113 + 1X4 + *i + ^2 + $1 + $2 + f^i + ^^2 + mirrors and permuted]. (57) 



The variables entering these expressions can be determined from the same linear system as above, 
expressed in terms of the new variables. For instance, the first equation becomes 




(58) 



mw 



The number of equations is the same, but we have fewer variables. The new sets {[1^2 x A/3]}, 
{[A/3 X N2]] and {As x As} have respectively 15, 15 and 25 variables, to which 20 extra variables are 
added for the diagrams in Figure 3 and their mirror images. So there are 88 variables, constrained by 
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81 linear equations, of which 73 only are linearly independent. The system is still underdetermined but 
unexpectedly allows to determine enough variables to compute the probabilities without further work. 

First one can show, by suitably combining the independent equations, that [A + B + C + C]/N is 
equal to a combination of local diagrams which turns out to be subdominant, of order for both an 
open or a closed boundary (order at the critical point, where m is the distance between iq and jo). 
As the four quantities are positive by construction, it means that each of them is at least of order i^, and 
can be neglected. Thus one relation determines four variables. Once these four variables are eliminated, 
one is left with a system of 72 independent equations for 84 variables. 

Being independent, the 72 equations allow to determine 72 combinations of variables. The point is 
that one can choose these 72 combinations in such a way that the probabilities can be fully expressed in 
terms of them only, thereby making the knowledge of the other 12 combinations useless. Alternatively, 
one may choose to solve the linear system for 72 variables, which then b(x;ome functions of the remain- 
ing 12. When inserted in the probabilities, all dependences in the 12 unknowns drop out completely. 
The set of the 12 variables that the system cannot determine is not unique, but a possible choice is 

-J-P, hP, W, OP + mirrors}. 

The counting of variables and equations is different when there are no insertions at other places than 
if) and jo since the mirrored equations are redundant. One finds that the linear system is again not 
invertible, but is nonetheless sufficient to compute all 2-site boundary probabilities. They have been 
computed in the massive model to the dominant order , which yields the universal terms. 

For an open boundary, we found that none of the diagrams in Figure 3 contribute to the dominant 
order, being at least of order t^. The probabilities P22,P23 and P33 have the same form t^[KQ{m^/t) — 
Kf){m^/t)\^ at dominant order, and only differ by their normalizations. These have been checked to be 
in agreement with the identifications obtained in Section 4. 

The case of a closed boundary is a bit more complicated. In this case the diagrams of Figure 3 
contribute to order (as we have seen above, none of the diagrams of Figure 2 contribute, irrespective 
of the boundary condition), and the probabilities read 

P22 = t^[ - ^KoimVtf - - ^) Kl^imVif + ^KoimVt)K(;imVi) 

-i^-B + J^) K^limVtr + K(mVt)K{mVt)) + . . . (59) 

P2S = e[-^K,{mVtf + (-^ - i^, + ^) K',{mVtf - ~ ^) K,{m^t)K'^{m^t) 

+ (|i - ^ + 4^) K'^irnVtf -{^-^) KimVmi^'imVi)) +... (60) 

^^33 = - ji^KoimVt)^ + ^K'o{mVif - ^Ko{mVt)K'^{msft) 

-{^ + ^) K^imVir + ^K'^{mVt)K'^'{mVi)) + . . . (61) 

Again they are in full agreement with the fields found in Section 4. 

We have also computed a few 3-site probabilities, when one of the insertion is a height 1 or a 
supercritical height. Then the same system as above can be used, the only difference is that the Laplacian 
has to be decorated by a local defect matrix, and only affects the calculation of the local diagrams. We 
have found for instance the connected probability P212 to have a height 1 and two heights 2 on a closed 
boundary, all separated by large distances, at the critical point (the expressions for off-critical 3-site 
probabilities are too long), 

'3 2\ / 1 1 1 (Gtt - 24)2 



-f- 



4 TtJ \mi2rnf3TO23 m^2^ti^'ii ™12™13™23 



67r-24 
+ 



1111 
+ 



mi2m^^m2^ m'(2m^^m2z mfj^iamlg m;p2'TT'i3TO23 

,-th J -th 



+ ... (62) 



where ruij is the distance between the i and j site. The connected probability is equal to P2i2,( 
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P212 — P2P12 — P22P1 — P21P2 + 2P|Pi. The previous formula for P212 is equivalent to that found by 
Jeng but allows for a more direct comparison with the field theoretic result, as the various terms 
correspond to specific Wick contractions. 
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